Statistics of Group Delays in 
Multimode Fiber with Strong Mode Coupling 

Keang-Po Ho, Senior Member, IEEE, Joseph M. Kahn, Fellow, IEEE 



(N : 

O ' 

<N : 

?H ■ 

Oh- 
< 



• I— I ■ 

> , 

Ph. 
q ; 

o . 

Oh' 



> 

o 



X 



Abstract — The modal group delays (GDs) are a key property 
governing the dispersion of signals propagating in a multimode 
liber (MMF). A MMF is in the strong-coupling regime when the 
total length of the MMF is much greater than the correlation 
length over which local principal modes can be considered 
constant. In this regime, the GDs can be described as the 
eigenvalues of zero-trace Gaussian unitary ensemble, and the 
probability density function (p.d.f.) of the GDs is the eigenvalue 
distribution of the ensemble. For fibers with two to seven modes, 
the marginal p.d.f. of the GDs is derived analytically. For libers 
with a large number of modes, this p.d.f. is shown to approach a 
semicircle distribution. In the strong-coupling regime, the delay 
spread is proportional to the square root of the number of 
independent sections, or the square root of the overall fiber 
length. 



I. Introduction 

MULTIMODE fiber (MMF) is widely used in short- 
range optical links |[l]-||3l, where it is often favored 
over single-mode fiber (SMF) because of relaxed connector 
alignment tolerances and reduced transceiver component costs. 
MMF supports propagation of multiple spatial modes having 
different group velocities, and thus different group delays 
(GDs), an effect called modal dispersion ID, ||5|. Even if 
a signal is launched into one spatial mode, bends, index 
imperfections and other perturbations cause the signal to 
couple into multiple modes |]4l, Q-ISl, making the signal 
subject to modal dispersion. Modal dispersion limits current 
commercial MMF links to 10 Gb/s per fiber up to about 300 m 
long m, m, and next-generation 100 Gb/s Ethernet MMF 
systems use ten fibers per link ifTOl . Techniques to increase 
the bit rate per fiber are desired. 

SMF, which is free from modal dispersion, is the dominant 
medium for longer transmission distances. Emerging long-haul 
systems use dual-polarization quaternary phase-shift keying 
and coherent detection to achieve a spectral efficiency of 2 
bits/s/Hz ifTTll . lfT2l . Usage of higher-order modulation formats 
ifTsl can at least double the spectral efficiency, but further 
increases are expected to become increasingly difficult |[T4l . 
because of limits posed by optical amplifier noise and fiber 
nonlinearity ifTSl . lfT6 l. Techniques to further increase spectral 
efficiency are desired. 
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Mode-division multiplexing (MDM) in MMF [T3, CSl, a 
form of multi-input, multi-output transmission, is a potential 
means to increase transmission capacity in both short- and 
long-distance optical networks. Like multipath propagation 
in wireless systems, the plurality of modes in MMF was 
long viewed as a strictly negative, bandwidth-limiting effect 
requiring mitigation, but is now seen as creating additional 
degrees of freedom in which to transmit information [:19l- 
1221 . Modal dispersion in MMF typically leads to a larger 
GD spread than that caused by chromatic dispersion. This GD 
spread determines the required cyclic prefix length in MDM 
systems using orthogonal frequency-division multiplexing 1211 
or the required number of equalizer taps in MDM systems 
using single-carrier modulation |22|. In other words, receiver 
complexity increases in proportion to the GD spread caused 
by modal dispersion. 

Effective mitigation of modal dispersion or optimal use of 
MDM requires a detailed understanding of modal dispersion, 
especially the effect of mode coupling on the modal GDs. 
Models for mode coupling were developed more than 30 
years ago ID, Q, JT], when MMF links used spatially and 
temporally incoherent light-emitting diodes. Virtually all the 
models ignore phase effects, and consider only power coupling 
between modes. Power coupling models are able to qualita- 
tively explain some observations, such as the scaling of delay 
spread with fiber length. Delay spread scales linearly with fiber 
length in the weak-coupling regime (e.g., short glass MMF), 
and with the square root of fiber length in the strong-coupling 
regime (e.g., plastic MMF) (23|. However, most modem MMF 
systems use spatially and temporally coherent laser sources, 
and power coupling models are not able to explain certain 
observations, such as a sensitivity of the impulse response to 
launched polarization ll24l . 

SMF supports propagation in two polarizations, and 
polarization-mode dispersion (PMD) has long been modeled 
using electric field coupling models ll25l . ll26l . Field coupling 
models have been used to demonstrate the existence of prin- 
cipal states of polarization (PSPs), which have well-defined 
group delays to first order in frequency Il26l - ll28l . In long 
SMFs, polarization modes are strongly coupled. In this regime, 
the differential group delay (DGD) between the PSPs scales 
with the square root of fiber length and follows a Maxwellian 
distribution Il27l - ll30l . PSPs form the basis of techniques for 
optical PMD compensation in direct detection systems. 

Recently, field coupling models have been extended to MMF 
131)1 . l32l . These models explain the polarization sensitivity 
of mode coupling and demonstrate the existence of principal 
modes (PMs). The PMs have well-defined GDs to first or- 
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der in frequency, and form the basis for optical techniques 
to compensate modal dispersion [[33l . The GD differences 
between PMs scale linearly with fiber length in the weak- 
coupling regime, and with the square root of fiber length in the 
strong-coupling regime ll32l . To date, however, the statistical 
properties of the GDs, which are of particular interest in the 
strong-coupling regime, have not been studied. 

Here, the statistics of the GDs are derived analytically for 
MMF in the strong-coupling regime, considering a number of 
modes ranging from two to infinitjQ. In this regime, regardless 
of the number of modes or the group delays in the absence of 
coupling, the GDs scale with the square root of fiber length or 
the square root of the number of independent fiber sections, 
similar to PMD in SMF in the strong-coupling regime ll27l . 
End-to-end modal dispersion effects are described, at each 
single frequency, by a random complex Gaussian Hermitian 
matrix or Gaussian unitary ensemble 12]. From such a model, 
the joint probability density function (p.d.f.) of GDs can be 
derived analytically. Here, closed-form expressions for the GD 
distributions are derived for small number of modes. For a 
large number of modes, the GD distribution asymptotically 
approaches a semicircle distribution with a radius or upper 
limit equal to twice its standard deviation. 

The remaining parts of this papers are organized as fol- 
lows. Sec. describes the random matrix model for MMF 
propagation. Sec. HHl provides closed-form analyses of the GD 
distribution in fibers with two to seven modes. Sec.lIVIpresents 
asymptotic expressions for the GD distribution in the limit of 
a large number of modes. Sees. |V] and |VT] are discussion and 
conclusions, respectively. 

n. Random Matrix Model for Multimode Fibers 

The propagation characteristics of a MMF, in particular, 
the local PMs and their GDs, can be considered invariant 
over a certain correlation length. Because bends, mechanical 
stresses and manufacturing tolerances induce mode coupling, 
the local PMs in sections separated by distance longer than 
the correlation length can be considered independent of each 
other. Throughout this paper, we consider the regime of strong 
mode coupling, where the total length of the MMF far exceeds 
the correlation length. The theory presented here is valid 
regardless of the actual correlation length or the modal GD 
profile within the correlation length. 

A MMF may be divided into K sections, with propagation 
in each section modeled as a random matrix. The length of 
each section should be at least slightly longer than the corre- 
lation length, so that the local PMs in the different sections 
can be considered independent. Although the approach used 
here is applicable even if each section has different properties, 
for convenience, we assume that all sections are statistically 
equivalent. 

In a MMF with D modes, modal propagation in the fcth 
section can be modeled as a x matrix M.'^^\ijj) as a 
function of frequency w. Here, we are only interested in the 

'Throughout this paper, "modes" include both polaiization and spatial 
degrees of freedom. For example, the two-mode case can describe the two 
polarization modes in SMF. 



Statistical properties of the modal GD, so for simplicity, we 
ignore any mode-dependent gain or loss. With strong mode 
coupling, the fcth section may be represented by the product 
of three D x D matrices 

= v('^'A"^■)(^^)^^*''^^fc- 1,...,^, (1) 

where ^ denotes Hermitian transpose, W''"' and V''"'-' are 
random unitary matrices representing the mode coupling at 
the input and output, respectively, and A''"''(aj) is a diagonal 
matrix describing the uncoupled modal GDs, i.e.. 



A('=)(w) = diag \e-^ 



(fc) 



(2) 



(k) 

where ' , k = 1, . . . , K, i = I, . . . , D, are the uncoupled 
GDs in the sections. 

In the absence of mode-dependent gain or loss, 
are all unitary matrices, such that 
M^^\uj)M^^''^ {uj) = I, where / is the identity matrix. With 
strong random coupling, both Z^'*^' and V''^-' can be assumed 
to be independent random unitary matrices, such that both 
input and output are randomly oriented. The model ([1} is 
similar to the matrix model of PMD described in ll30l . Il35l . 

The model here is valid regardless of whether or not the 



vectors t 



•(fc) 



' ~2 



'D 



1, . . . ,K, have 

•(fe) 



the same statistical properties. The vector f'- ' may even 
be a deterministic vector, identical for each section. For 
convenience and without loss of generality, we assume that 
J2i T^^^ ~ 0, i.e., we ignore the mode-averaged delay of each 
section, as it does not lead to modal dispersion. 

Using \sk) and \tk) to denote the input and output modes, 
respectively, we have 

=7WWH|sfc) (3) 

and \sk) = M'^^^\uj)\tk) ■ Similar to the analysis of PMD 
EH, ESI, ll35l . the GDs correspond to the eigenvalues of 
jM\l;^M'^^^'^ where M^^^ = dM'^^\uj)/duj. With only a 
single section, we may verify that 

jM'^^^M^^^'^ ^V'^^^T^^^V'-^^'^ (4) 

with V'^'^-' as the local PMs in the fcth section, and where 



T'^^^ = diag 



^k) (k) 



... J^) 



(5) 



is a diagonal matrix of their GDs in the fcth section. With 



0, we have tr (T^''^) = and 

= 0. 



tr 



(6) 



Physically, the ith local PM experiences an uncoupled GD 
Tj''"^ without mixing with other modes. Because the diagonal 
matrices T*-'^' are real matrices, all matrices jA^if 'AI''^-'^, 
fc = 1, iiT, are Hermitian. 

When K sections of MMF are cascaded together, the overall 
propagation matrix becomes 

TWW =7WW7W(^-i)...7W(2)x(i). (7) 

The overall PMs and their GDs correspond to the eigenvectors 
and eigenvalues of ll3T1 . Il32l 



(8) 
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Because 

+ (9) 

we obtain 

• • • + Wx(^-i) . . . 7W(2)_y^(i)7v^(t)t^ (^Q) 

From ( [Tol l, the overall matrix j.mL*''7M'-*-'^ is the sum- 
mation of K random matrices. All those K matrices have 
eigenvalues identical to those of ^ and (|5]i. However, their 
eigenvectors are independent of each other The first matrix 
jA^o/^^A^'^''^ has eigenvectors derived from V'^^K The sec- 
ond matrix jM'^^'>M^~^^M'^^~^'>''M^^'>'' has eigenvectors 
derived from 7W('^)v('^"i). Matrices V^^^ and A^(A')v(^"1) 
are both unitary matrices and are obviously independent of 
each other All the K matrices summed to form Q are inde- 
pendent of each other with eigenvalues given by the vectors 
T^''\ k ~ l,...,K. Even for the case that all vectors t''"'-' 
are deterministic and identical, all the K matrices summed to 
form Q are independent, owing to the different directions of 
their independent eigenvectors. 

The matrix elements of Q should be identically distributed 
Gaussian random variables from the central limit theorem 
(CLT). The matrix elements of Q, gij, i,j ~ 1,...,D, are 
the summation of K identically distributed random variables, 
as seen from ( fTOl i. If K is very large, gij are Gaussian random 
variables from the CLT. Because all K component matrices in 
([Tol l are Hermitian, ^? is a Hermitian matrix. The diagonal 
elements g^^, i = 1,...,D, are all real Gaussian random 
variables with variance equal to cr^. All non-diagonal elements 
9i,j, i 7^ are complex Gaussian random variables with 
independent real and imaginary parts, which have variance 
equal to o-g/S. Thus, the elements g^j-, i ^ j has variance 
cTg. The value cr^ depends on the number of modes D, the 
number of sections K, and the variances of the uncoupled 
GDs described by t^^\ If the D vectors in V'^'^^ given by (|4]i 
are assumed to be independent of each other, it can be shown 
that 

k=l k=l 

where cr^(fc) ,k = l,...,K, are the variances of the GDs in the 
sections. If all K sections have the same modal GD profiles, 
we have 



where cr^ are the GD variances in all sections. 

In random matrix theory, the matrix Q is described as a 
Gaussian unitary ensemble |J2], Sec. 2.5]. Typically, a Gaussian 
unitary ensemble does not have any constraint aside from the 
variance of its Gaussian elements. However, in (fTOl i. the matrix 
components have zero trace so that 

tr(g) = 0. (13) 



In other words, is a zero-trace Gaussian unitary ensemble. 
The GDs in a MMF are statistically described by the eigen- 
values of the zero-trace Gaussian unitary ensemble. 

The assumption that the diagonal and off-diagonal elements 
of Q have the same variance (fTSl i is valid only if the 
D orthogonal vectors in V'^'^' are independent of each other 
However, the condition of orthogonality implies that the Dih 
vector is determined by the other D — 1 vectors. Using 
numerical simulation, we have found that all diagonal elements 
of Q have equal variance of (1 — D^^)a^, all off-diagonal 
elements of Q have equal variance of [1 + 1/D{D — l)]cr^, 
and the average variance of all elements of Q is ( fT2l ). 
The theory of ||2] Sec. 14.3] is able to describe a random 
Hermitian matrix in which different elements have different 
variances. Numerical simulations of the zero-trace Gaussian 
unitary ensemble with unequal variances that is considered 
here show no observable differences from analytical results 
derived assuming all matrix elements have equal variance ||2] 
Sec. 3.3]. 

III. Modal Dispersion in Few-Mode Fibers 

In the regime of strong mode coupling, the PMs and their 
GDs are given by the eigenvectors and eigenvalues of the zero- 
trace Gaussian unitary ensemble described by (ITOt and ([13). 
Without loss of generality, after normalization, the elements 
of Q may be assumed to be zero-mean identically distributed 
Gaussian random variables with variance cTg — 1/2, similar to 
the classic normalization of Mehta ||2]F1 Before normalization, 
(Tg is given by either ( fTTT i or (fT2t . The diagonal elements of 
Q are real with a variance of cr^ ~ 1/2. The off-diagonal 
elements of Q are complex Gaussian distributed with inde- 
pendent real and imaginary parts, each having a variance of 
1/4. From (fT2] |. each section can be taken to have a normalized 
standard deviation of GD given by (t^ = \/D/2K. With this 
normalization, the notation in this section is similar to that in 
Mehta El. 

A. Joint Probability Density 

The joint p.d.f. for a Gaussian unitary ensemble without the 
zero-trace constraint is well-known. The ordered joint p.d.f. of 
the eigenvalues a Dx D Gaussian unitary ensemble is ||2] Sec. 

3.3] na 

ao n (A.-A,fcxp(-^An , (14) 

D>i>j>0 \ 1=1 / 

where the eigenvalues possess the order constraint Ai < A2 < 
■ • • < and are all real valued, and an is a constant such 
that the joint p.d.f. integrates to unity. The eigenvalues are 
the normalized GDs. An analytical expression for ajj can be 
found in |2j Theorem 3.3.1]. Because a permutation with a 
different ordering of eigenvalues is equivalent to any other 
permutation, the unordered joint p.d.f. is just 1/Dl of ( IS.lt 
but without the order constraint [2, ch. 5]. 

^In some mathematical literature, which is relevant to Gaussian orthogonal 
ensembles but not Gaussian unitary ensembles, the matrix elements are 
assumed to have unit variance. Proportionality constants are ignored in some 
of those references. 
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With zero trace, 

trig) = Xi + X2 + --- + Xd ^0, 
the ordered joint p.d.f. of Ai, . . . , Xd-i becomes 

Pd{Xi,...,Xd-i) = Pd Yl (^i-^j)^cxp 

D>i>j>0 

with the order constraint 

Xi<X2<---<Xd, 



(15) 



(17) 



where the constant is determined by requiring (IS.2t to 
integrate to unity, but is not the same as the an in (IS. 11 1. 
The unordered joint p.d.f. is the same as (|S.2| i but is a factor 
of 1/Dl smaller and does not have the order constraint (ITTt . 
The statistical properties of GD are fully specified by the joint 
p.d.f. (|S?2] | with the constraint (flST l. 



B. Two-Mode Fiber 

Two-mode fiber is the simplest case, and may correspond 
to the two polarization modes in a SMF, i.e., the well- 
known PMD problem. The purpose here is not to derive new 
properties of PMD, but to verify that the general random 
matrix model is applicable to PMD. 

With A2 = -Ai, the p.d.f. (|S^ for £» = 2 becomes 



P2(Ai) =/324A2e 



(18) 



As in the PMD literature, we define Ai^2 = ^^'t ^nd find 
/32 — obtaining 



-r^e-^'/^ T > 0, 

TT 



(19) 



which is the well-known Maxwellian distribution with nor- 
malized mean DGD of f = 2\plpK = 1.60. Random matrix 
models specialized to the two-mode case were used to derive 
the Maxwellian distribution in 1281 . Il30l . The second moment 
of both Ai and A2 is 3/4. 



C. Three-Mode Fiber 

Three-mode fiber is the next-simplest case, with A3 = 
-(Ai + A2). The joint p.d.f. dO for Z) = 3 of ([iTli is 



P3(Ai, A2) =/?3(A2 - Ai)^(2Ai + A2)'(Ai + 2A2)' 

X e"^i"^2-(^i+-^2)'^ (20) 



Without the order constraint (|T7t . the marginal p.d.f. of the 
GDs is 



1 /■+°° 

P3a(A) = - / P3(A, A2)dA: 
6 J -00 



(21) 
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Fig. 1. Statistics of tlie nomialized GDs for thiee-mode fiber, including 
the marginal p.d.f. P3a(^). the p.d.f. ps\-^ (Ai) of the smallest delay, and the 
p.d.f. P3A, (A2) of the middle delay. To facilitate comparison, both p^Xi (-^i) 
and 'P:iX;{^2) are scaled by a factor of 1/3. 



Some algebr^ yields the constant jS^ = Ay/S/ir and 



P3a(A) = 



^y|(27A4-18A2 + 5)cxp(^-^A2). (22) 

Figure [T] plots the marginal p.d.f. P3\{X), which exhibits 
three peaks, corresponding to the values where Ai, A2, A3 are 
concentrated. The p.d.f. P3a(A) is symmetrical with respect to 
A = due to the symmetric nature of the three eigenvalues 
for the 3x3 random matrix Q. The middle eigenvalue A2 is 
concentrated near zero. The variance of A is cr^ = 4/3. 

The p.d.f. of the eigenvalue Ai, corresponding to the 
smallest delay, can also be found using ( IS.2b with the order 
constraint (fTTT i. The condition Ai < is required in order to 
conform to the zero-trace constraint (ITSt . The p.d.f. of the 
smallest eigenvalue is 



P3Ai(Ai) 



P3(Ai, A2)dA2, 



(23) 



Ai 



P3Ai(Ai) = 3p3A(Ai)erf 



3|Ai| 
^/2 

9V3Ai(3A? - 5) 



Stt 



■exp( 



-6A?), 
Ai < 0. 



(24) 



Due to symmetric nature of Ai and A3, we have p^Xs (A3) = 
PsAi (^A3). Figure[T]also shows ^P3Xi (Ai) where the p.d.f. is 
scaled by a factor 1/3 such that P3Ai(Ai) given by ( |24] | is 

^ All the calculations perfoiTned for D > 3 require three steps. In the first 
step, A|+A| + - ■ ■+A|,_-^-|-(Ai+A2 + - • ■ + A_D_i)'^ is linearly transformed 
to D/{D — 1) X A|+x| + - ■ In the second step, the linear transfomi 

is substituted into YlD>i>j>o(^^ ~ "^j)^ ■ which may be expanded to a 
summation of terms in the foiTn of c/j./,,..,/^_j A{^a;2^ ■ • x^^^^^, where 
fi are the exponents and c/i,/2,.../£j_i the con'esponding coefficients. 
The last step is the integration over X2, x-^, . . . , X]j_,i using |37| Sec. 
3.461]. These calculations ai'e tedious, but can be performed using symboUc 
mathematical softwai'e, such as Maple or MuPAD. 
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TABLE I 

Statistical parameters for the normalized GDs of a 
three-mode fiber 



p.d.f. 


Mean 




Variance 






m(A) 







4 
3 






P3Ai(Al) 


27 
16 


/T 


4 1 9v^ 
3 ^ Sir 


729 _ 
1287r 


0.1407 


P3A2(->*2) 







4 9^ 
3 4ir 


= 0.0928 






27 / 
16 V 


T 

n 


0.1407 






P3(A3-Ai)(A) 


27 / 
8 V 


Y 

IT 


4_ 729 _ 
327r 


1 27^3 _ 
47r ~ 


0.4700 



nearly the same as P3\{X) given by ( |22] ) near the first peak of 
P3a(A). 

Similarly, the p.d.f. of the middle eigenvalue A2 can be 
found by 



P3A2(A2) 



-\\2\ 



Psi^i, A2)dAi, 



(25) 



or 



P3A2(A2) = 3p3A(A2)erfc 



3IA2I 



V2 

9V3|A2|(3Ai-5) 



cxp(-6A2), (26) 



Comparing (|22] |. ( l24l i. and (|26] |. the marginal p.d.f. P3a(A) 
is found to be the combination of P3Ai(Ai), P3A2(A2), 
and p3A3(A3): 



P3a(A) = - b3Ai(A) +P3A2(A) +P3A3(A)] 



(27) 



Figure [U also shows ■|p3Ai(Ai) and ■|p3A2(A2) which are 
concentrated near the corresponding peaks of the marginal 
p.d.f. P3a(A), confirming that each peak of the marginal 
p.d.f. P3a(A) corresponds to an individual eigenvalue. The 
statistical parameters of the normalized GDs of a three-mode 
fiber are presented in Table |I] In Fig. [T] the peak for the 
middle delay A2 is narrower than those for the maximum and 
minimum delays Ai and A3. In Table |l] the variance of A2 is 
smaller than the variance of Ai or A3. 

The difference between the maximum and minimum eigen- 
values is the normalized delay spread of the MMF The 
p.d.f. of the delay spread for three-mode MMF is 



P3(A3-Ai)(A) 



A/3 



P3 



-A/3 



A- A, 



A2 dXo 



(28) 



or 



P3|A3-A.)(A) = |2(A» -9A')cxp (^-jA'j 

+ iyf(x«-6A. + m')„p(-lA')o,f(^ 

A > 0. (29) 

The statistical parameters of the delay-spread are also given 
in Table U For three-mode fiber, many properties of the 
eigenvalues or normalized GDs can be computed analytically 
in closed form. 




-202 
Normalized Group Delay, A, 

Fig. 2. Tlie marginal p.d.f. of the normalized GDs for a four-mode fiber, 
comparing analysis and simulation. 



D. Four-Mode Fiber 

Four-mode fiber is a particularly simple case beyond the 
two-mode fiber, as it represents a fiber with two spatial modes 
and two polarizations. Fibers with two spatial modes have 
been used for dispersion compensation ||38l , ||39l and for 
transmission experiments 11211 . POl . We note, however, that a 
weakly guiding fiber with circular core cannot support exactly 
two spatial modes BTI . Although the relatively short distances 
used in 11211 . 1401 may not be sufficient to ensure strong mode 
coupling, the strong-coupling regime can be expected in future 
long-distance transmission systems. 

The p.d.f. ( IS.2I 1 for D ~ A becomes 

P4(Ai,A2,A3)-/34 n (A«-Aj)2e-^'-^^-^3-A^, (30) 

with the zero trace constraint A4 = — Ai — A2 — A3. 

Using the unordered joint p.d.f., the marginal p.d.f. of A is 



P4a(A) = 



+00 p-\-00 



P4(A,A2,A3)dA2dA3. (31) 



After some calculations, we obtain 

2v^e-4/3AV4096,fi 1024 



729 



^ +81^ 



_5_ 

81 



^^"^^ ^ \/tt V6561 i 01 o 

02) 

which has a variance of 15/8. 

Figure |2] shows the marginal p.d.f. of the normalized GDs 
in a four-mode fiber given by (l32t . The marginal p.d.f. has 
four peaks, corresponding to the GD of four different PMs. 

To verify the marginal p.d.f. in Fig. |2] the modal dispersion 
of a four-mode MMF has been simulated. The fiber has 
K = 256 independent sections. In each section, the four 
modes are chosen to have deterministic delays of +t, +t, — r, 
and — r where r = K to ensure that the elements of Q 
have a variance of = 1/2. This particular choice could 
describe a fiber where, in each section, the DGD between 
the two polarization modes is negligible compared to that 
between the two spatial modes. The random unitary matrices 
and V^*-'), k = l,...,K, are first initialized by 4 x 4 
random complex Gaussian matrices and then converted to 
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unitary matrices using the Gram-Schmidt process B2] Sec. 
5.2.8]. All sections have independent matrices U'^'^^ and V'''"-'. 
A total of 400, 000 eigenvalues are used in the curve shown 
in Fig. 2. 

In Fig. |2l the simulation results show excellent agreement 
with the analytical p.d.f. ( |32] i. Although the modes in each 
section have only two GDs, with strong mode coupling, a 
p.d.f. having four peaks is obtained. In the strong-coupling 
regime, similar results would be obtained using any uncoupled 
GDs in each section, provided that the four GDs sum to zero 
and have a variance of 2/K. For example, the four modes 
may have modal delays of +ti, +ti, +ti, and — 3ti with ti = 
^/IjiK. Similar results may also be obtained if the GDs in 
each section are, for example, +r2, +T2, — T2, — T2, where 
follows a statistical distribution with second moment 21 K. 

The marginal p.d.f. of the smallest or largest eigenvalues Ai 
and A4 may be found by suitable integration of dSOl l. Unlike the 
case of three-mode fiber, it does not seem possible to obtain 
closed-form expressions for the individual marginal p.d.f.'s of 
the ordered eigenvalues. 

E. Other Few-Mode Fibers 

The marginal p.d.f.'s of the GDs in MMF with larger 
number of modes may also be obtained analytically. Following 
the above procedure, the marginal p.d.f.'s of fiber with five, 
six, and seven modes are 



P5a(A) = 



^g-5/4A^ / 78125 
196608 



15625 
8192 



24375 
' 8192 



A^- 



1975 
2048 



903 
4096 



, (33) 



P6a(A) 



30e-6/5A / 13436928 10 4478976 8 

A — „ A 



V 244140625 



9765625 



2581632^, _ 102816^, ^m^^,^^^^ (34) 



1953125 



78125 



15625 15625 J 



and 



^426-^/^5^" / 1977326743 
= ^ I ^TT^r^rT^rT^T^ A 



282475249 
1632586752 
17707375 



A^'^ + 



V 146932807680 

98001617 8 47883143 g 
A — A 



120932352 
212219 



30233088 
88175 \ 



13436928 746496 1492992 / ' 



(35) 



respectively. The variances of these distributions are 12/5, 
35/12, and 24/7, respectively. For all cases of D from two 
to seven modes, the variances are given by ^{D — D^^), a 
reduction by a factor of 1 — compared with the case 
without the zero-trace constraint. 

Figure [3] shows the marginal p.d.f. of the normalized GDs of 
fibers with five, six, and seven modes. The number of peaks 
in the marginal p.d.f. is the same as the number of modes. 
In general, the peaks closer to A = are both higher and 
narrower than those farther from the origin. Those peaks cause 
ripples to appear in the marginal p.d.f., and the ripples are 
still significant in seven-mode MMF. As the number of modes 



0.25 - 



■q 0.15- 




-5 -4 -3 -2 -1 1 2 3 4 5 
Normalized Group Delay, ^ 

Fig. 3. The marginal p.d.f. of nonnalized GDs for fibers with five, six, and 
seven modes. 



increases, the peaks in the marginal p.d.f. move closer together 
and merge. As the peaks become indistinguishable for fibers 
having many modes, the marginal p.d.f. should approach a 
limiting distribution. 

Numerical simulations, similar to those in Fig. |2] have been 
used to verify the analytical p.d.f.'s of Fig. |3] In all the cases, 
simulation and theory match with each other 

Numerical simulations have been conducted to further verify 
the variance reduction factor of 1 — D~^, which is most 
significant for D < 3. Random realizations of zero-trace 
matrices Q of the form ( fTOl i exhibit no observable variance 
reduction, although the empirically estimated p.d.f.'s, when 
scaled by this reduction factor, are found to match with ([19) 
and ( |22] |. 

The reduction of variance may be seen as related to de- 
grees of freedom. A random Hermitian matrix without zero- 
trace constraint has degrees of freedom, corresponding to 
D{D—l)/2 complex off-diagonal elements and D real diago- 
nal elements. The zero-trace constraint reduces the degrees of 
freedom by one, proportionally affecting a fraction 1/D^ of 
matrix elements. The zero-trace constraint from (IS.l) to (IS.2l i 
reduces D degrees of freedom to 13 — 1 degrees of freedom, 
proportionally affecting a fraction 1/D of matrix elements. 
While analytical results scaled to the same variance are 
consistent with numerical simulations, the variance reduction 
factor of 1 — D^^ requires further study. 

In the simplest case of D = 2, a zero-trace Gaussian unitary 
ensemble can be generated numerically by three methods. The 
first method is based on Q given by the summation ( fTOb . for 
example, with <t^ = 1/2. The second method is based on a 
random 2x2 Hermitian matrix A but with 02,2 replaced by 
— ai 1. The third method is based on generating random 2x2 
Hermitian matrix A, finding its eigenvalues Ai and A2, and 
selecting those with |Ai + A2I smaller than a small number 
In the second and third methods, the elements of the matrix A 
have variance of 1/2. The first two methods give eigenvalues 
with the same variance but the third method gives eigenvalues 
with a variance 3/4 time smaller than the first two methods. 
In the third method, the variance of the diagonal elements 
are 1/4 and the off-diagonal elements have variance of 1/2. 
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The zero-trace constraint on the eigenvalues (fT5t reduces the 
variance of the matrix elements but the p.d.f. maintains the 
same shape. 

The variance reduction of the diagonal elements is due to the 
zero-trace constraint ( fTTt . which selects those matrices with 
smaller diagonal elements (in the general case, a factor of 
1 — D^^ smaller). The average variance among all elements 
is a factor 1 — D^^ smaller, the same as the reduction factor 
for the variance of the eigenvalues. 

As the average variance for all matrix elements of Q is o-^ 
given by ( fT2b . for all MMF studies in this section, and without 
normalization, the variance of the GD is 

4d = K<yl (36) 

with (T^ defined by ( fT2b . If the GD is characterized by its 
standard deviation <Tgd, it is always proportional to the square 
root of the number of independent MMF sections. 

IV. Modal Dispersion in Many-Mode Fibers 

With a large number of modes, a Gaussian unitary ensemble 
without the zero-trace constraint is described by a semicircle 
distribution with radius \/2D ||2] Sec. 4.2]. With the normal- 
ization used in Sec. [nil the variance of the eigenvalues is 
D/2. This semicircle law was first derived by Wigner for 
large random matrices P3l . ll44l . The Wigner semicircle law 
is universally valid for many different types of large random 
matrices B31 . 1461 . A Gaussian unitary ensemble, even with 
the zero-trace constraint (fTTt . should follow the semicircle 
distribution. As an alternative to considering as a Gaussian 
unitary ensemble, a more straightforward derivation using the 
CLT for free random variables is given in a later part of this 
section. 

In free probability theory, free random variables are equiv- 
alent to statistically independent large random matrices ||47l , 
PHI . The CLT for the summation of free random variables 
gives the semicircle distribution BTI . P9l . The matrix Q (fTOl i 
is the summation of many independent random matrices. The 
CLT for free random variables states the following: Let X^, 
k = 1, . . . ,K, be identically distributed independent zero- 
mean free random variables with unit variance. The summation 



yii = 



Xi + X2 + ■ ■ ■ + Xj 



K 



K 



(37) 



is described by semicircle distribution with radius of two and 
unit variance 



kl < 2 



'0 otherwise 



(38) 



as K approaches infinity. 

In the above theorem, when free random variables are 
represented by large random matrices, the distribution of a 
free random variable is equivalent to the distribution for the 
eigenvalues of the random matrices. When the CLT of free 
random variables is applied to Q given by ( fTOl i. if the variance 
of the zero-mean GD per section is cr^ for all K sections, 
the eigenvalues of Q are described by a semicircle distribution 
with radius 2\lKa^ and variance Ka^. Equivalently, the GD 
of the MMF has a semicircle distribution with variance Kai. 




-1 1 

Normalized GrouD Delay, "k 

(a) L> = 16 




-1 1 

Normalized Group Delay, ^ 

(b) Z) = 64 




-1 1 

Normalized Group Delay, ^ 

(c) D = 512 

Fig. 4. Simulated marginal p.d.f. of nonnalized GDs (ripply blue curves) 
compared with semicircle distribution (smooth red curves) in MMF with (a) 
16, (b) 64, and (c) 512 modes. 



Note that the normalization used in this section based on 
the eigenvalues of Xu and yK in ( |37] | is customary in free 
probability theory. However, the normalization used in Sec. Hill 
is based on the matrix elements of Q, similar to that in Mehta 

El. 

Figure |4] compares the simulated marginal p.d.f. of MMFs 
having D = 16, 64, and 512 modes to the semicircle 
distribution. Each MMF is comprised of K = 256 sections. 
In each section, the GDs are deterministic, with the first Z5/2 
modes with a delay of t and the other D /2 modes with a delay 
of — r. The simulated curves are obtained from 1, 600, 000 and 
640, 000 and 102, 400 eigenvalues for D = \% and 64 and 512, 
respectively, using a step size of 0.025 along the A axis. The 
normalization t = is made to facilitate comparison 

with a semicircle distribution with radius of 2. The model 
here is valid as long as the component matrices in (fTol i may 
be modeled as free random variables IBtI . 

In Fig. m the simulated distributions match the semicircle 
distribution well for D = 64 and 512 modes. For a fiber 
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having D = 16 modes, the distribution is close to a semicircle 
distribution, but has an obvious periodic structure with 16 
peaks. The ripples become less obvious as D increases from 
16 to 64 to 512. Upon close examination of the curve for 
D = 64, the ripples seem periodic, similar to those in MMF 
having D = 16 modes, but much smaller. 

The semicircle distribution, describing the GDs in a MMF 
with an infinite number of modes, has strict upper and lower 
limits, and thus a strictly bounded GD spread. In designing 
systems for MMF with a finite but large number of modes D, it 
will be sufficient to provide a GD tolerance just slightly larger 
than the maximum GD spread of the semicircle distribution, 
which is given by 4crgd. 

The GD relationship ( [36] l remains valid when the number 
of modes D is vary large. With a large number of modes, 
the relationship a^^ = Ka^ can be derived directly from free 
probability theory. 

V. Discussion 

The scaling of modal dispersion with fiber length in MMF is 
similar to the scaling of PMD in SMF In a MMF shorter than 
the correlation length over which the local PMs can be consid- 
ered constant, the GD increases linearly with fiber length. By 
contrast, in a MMF much longer than the correlation length, 
the number of independent sections K is large, and strong 
mode coupling can be assumed. In the strong-coupling regime, 
a parameter describing GD per unit length may be defined as 
fkm = cTr/^/L^, where Ls is the fiber length per section, 
measured in kilometers. The overall GD, if characterized by 
(Tgd given by ( |36] |, is equal to ^/L^ay^ni, where Lt = KLs is 
the total fiber length. 

In practice, there are advantages to introducing strong mode 
coupling in order to reduce the modal delay spread. In direct- 
detection systems, this can reduce intersymbol interference, 
whereas in systems using coherent detection, this can reduce 
the temporal memory required in digital compensation of 
modal dispersion. Recent MDM experiments 1211 . 1221 . l40l . 
performed in short spans of MMF, were probably not in 
the strong-coupling regime. Future long-distance systems are 
likely to be operated in the strong-coupling regime, especially 
if strong mode coupling is used to reduce the overall GD 
spread. In MMF, spatial mode coupling is governed, in part, 
by mode groups |50j |. Typically, coupling between modes in 
different groups is weak, with coupling length as long as 25 
km |l6l, while coupling between modes in the same group is 
strong, with coupling length less than 1 km ITtI . In order 
to reduce the GD spread in MMF, coupling between modes 
in different groups should be enhanced. In manufacturing of 
SMF, spinning is used to reduce the polarization coupling 
length below 100 m, thereby reducing the DGD due to PMD 
II25I . Manufacturing processes for MMF may perhaps be 
modified to increase spatial mode coupling in order to reduce 
the GD spread. 

As seen in Figs. [3] and St a), in the marginal p.d.f. of GD, 
the number of peaks is the same as the number of modes, 
and the separation betweens adjacent peaks (relative to the 
semicircle radius) decreases with an increasing of number of 



modes. In the absence of the zero-trace constraint (ITSt . ripples 
can be observed observed in Gaussian unitary ensembles up to 
at least 51 x 51 ||2l Fig. 6.1]. With a zero-trace constraint, the 
ripples are larger than those without the constraint. In Fig. HJb) 
with D = 6A modes, ripples are observable and seems to be 
very regular. As the number of modes increases, the ripples 
becomes narrower, similar to the Gibbs phenomenon llSTI . Il52l 
for Fourier series. 

Higher-order modal dispersion effects are outside the scope 
of this paper. In higher-order modal dispersion, the PMs and 
their GDs can vary with frequency ||53l. These effects are 
analogous to polarization-dependent chromatic dispersion and 
depolarization observed in SMF with PMD EH, ||55]. In the 
case of SMF with PMD, the properties of PMD to arbitrary 
order depend on a single parameter In the case of MMF 
with modal dispersion and strong coupling, the higher-order 
properties of modal dispersion depend only on the number of 
modes and a single parameter, which may be taken to be the 
GD standard deviation (Tgd given by (l36l i. 

In this paper, we have studied the distribution of GDs, but 
not the impulse response of a MMF. At a single frequency, 
the impulse response of a D-mode fiber consists of D narrow 
pulses with GDs described by the distribution (IS.2l l |!32l, and 
with weights depending on the PMs excited by the transmitter 
launch conditions. Considering a modulated signal occupying 
a finite bandwidth, because of higher-order effects [531, those 
D narrow pulses broaden and may merge with each other The 
overall duration of the impulse response is described by the 
duration of the p.d.f. of the GD, as shown Figs. [T] to |4] In 
a fiber with many modes, where the p.d.f. is the semicircle 
distribution shown in Fig. |4] the impulse response duration is 

just A(7gd. 

VI. Conclusion 

In the regime of strong mode coupling, a MMF may be 
modeled as the cascade of many independent sections, which 
are described by statistically independent random matrices. 
The GDs are given by the eigenvalues of Gaussian unitary 
ensemble with zero-trace constraint. Marginal p.d.f.'s of the 
GDs in fibers with two to seven modes have been derived an- 
alytically. Numerical simulations of the p.d.f.'s are in excellent 
agreement with analytical results. In a fiber with many modes, 
the GD is shown to follow a semicircle distribution from free 
probability theory. Numerical simulations have been conducted 
for fibers having D = 16, 64, and 512 modes to compare to 
the semicircle distribution. 
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Clarification and Supplement 

Keang-Po Ho 

The main purpose of this note is to clarify the explanations in the paper [T]. The group delay statistics is also extended to 
multimode fiber with more than D ^ 7 modes, the limitation of the direct method in H]. 

For the random Gaussian unitary matrix Q without trace constraint, the joint probability density function (p.d.f.) of the 
eigenvalues is given by 12] Sec. 3.3] ||3l 



Pnc(x) 



■ exp 




(S.l) 



an Y\. 

D>i>j>0 

where an is a constant ID Eq. 3.3.10] such that the p.d.f. integrated to unity. Here, the notation of the eigenvalue changes 
to Xi instead of Ai in the paper, consistent with l2| and other literatures. The zero-trace random Gaussian unitary matrix is 
equivalently the statistics of Q — jyt^'G with joint p.d.f. 

Pzt(x) = I3dS 




n 

-D>i>j>0 



D 



(S.2) 



The constant Pu is not the same as an and can be found equal to (3d = V TrDan in later part of this note. 

For the p.d.f. Pnc(x) (IS. II ) without constraint, the elements of the random Gaussian unitary matrix all have the same variance. 
For the p.d.f. Pzt(x) (IS. 2b with zero-trace constraint, the diagonal elements are a factor 1 — 1/D smaller in variance than other 
off-diagonal elements. The combined variance of all elements of Q — jjt^Q is a factor 1 — 1/ smaller than that for Q without 
constraint. The eigenvalue variance of Pzt(x) iS.2i with zero-trace constraint is a factor 1 — smaller than that of Pnc(x) 

( IS.lt without constraint. The derivation in |[T] is consistent with this simple explanation but not as elegant here. 

The fiber simulation of HI, ID equivalently generates zero-trace Gaussian random matrix in which the variance of each 
element is not reduced by the factor 1 — l/D^. The random matrix in |[T] is equivalently 



D 



that the average variance in all elements is the same as Q without constraint. In [|4l, the zero-trace constraint is directly in 
Q without reducing the variance of the diagonal elements and the modification hy Q — §lrC7 is not required.. In [I], the 
p.d.f. derived from (IS.2l i is required to scaled up by the factor of D/V D"^ — 1 to match the simulation results. In ID, the 
eigenvalue p.d.f. is pre-scaled to unity variance. 



The eigenvalue distribution for zero-trace Gaussian unitary ensemble is derived according to 

/ D \ 



Pd{xi) = (3_ 



D 



/+00 /'+C30 / 

-oo J-oo \^^^ 



D>i>j>0 \ i=l 



dx2 ■ ■ ■ dx 



D- 



(S.3) 



First of all, D x D Vandemonde determinant gives 



det 



ij = l,2....-D 



r>>i>j>o 

where det[ ] denotes a determinant. Follow the method of ||2] Ch. 4] and directly from the properties of determinant, the 
Vandemonde determinant can be expressed by the Hermite polynomial as 



det 



i,j = l,2,...D 



= dot 



Hj^i{xi) 



i,j=l,2,...D 



and 



det 



dot 



i,j=l,2,...D 



^Hj-i{xi + Ci) 



i,j=l,2,...D 



(S.5) 



(S.6) 



where Hn{x) is the Hermite polynomial as given in [21 and q is constant. The leading terms in both 2 and 



are x;^ . Using the Hermite polynomial, the p.d.f. ( IS. Il l can be expressed as IS) 

Pnc(x) = det [KD{x^,Xj\^^-^ .^ JJ , 



where 



D-l 



k=0 



H^{x)Hk{y)e 



^x^/2-y-/2 



(S.7) 



(S.8) 



Without constraint, the eigenvalue distribution is given by j^Koix, x) and Kd{x,x) is called correlation function ||2]. 
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X 



Fig. 5. The distribution of Dpr){x) (solid curves) for zero-trace Gaussian unitary ensemble as compared with the correlation function of Kr){x,x) (dot 
curves) without constraint. 



With zero-trace constraint, the p.d.f. 

Pzt(x) 



becomes 

Pd 
Dar 



6 Xi^ dot [KD{xt,Xj)]^ j^-^^ 



2....D ■ 



Similar to ||5|, IS) but using Fourier instead of Laplace transform, 



Pzt(x) 



or 



Pzt(x) 



2Tr Da o j_ 

a r+ao 
PD 



r+oo f " \ 

J dot [A'D(a;,,a;j)]^_^.^^_2,...£)'^^P ( ^^^^^^^ 1 dw, 



2iTDaD J- 



det 



lUJ lUJ 
KD[Xi + —,X.j + — 



exp [ I dti^ 



(S.9) 



(S.IO) 



(S.ll) 



ij = l,2,...D 

changes from Xi to Xi + iuj/2 using ( IS.6I 



In ( IS.l II ). the argument for the Hermite polynomial in 

Similar to the integration of (IS.7t to obtain the eigenvalue distribution, the p.d.f. dS.SI l is given by 

f + 00 



Poix) 



Pi 



271 DaD 



IbJ ILjJ 

Kd [x + —,x+ — 



exp 



D 



or 



Pd{x) 



/3. 



D 



271 Dao 



+ 00 



exp 



D 



D-l 

E 

fe=0 



Hi X 



-duj. 



(S.12) 



(S.13) 



In the integration (IS. 13b . we may conduct the first integration over x before integration over uj to obtain Po/ctD = x/ttD. The 
integration pd{x) dS.Ob becomes 



Pd{x) 



1 



exp 



D 2 



X+-U 



Hi {x 



J k=0 



2^kUf^ 



-do;. 



(S.14) 



Using the substitute of s = iuj, the integration ( IS. 14b is very similar to the Mellin's inverse formula for Laplace transform 
with an integration from x — zoo to x + ioo. The following algebraic expression can be obtained based on Laplace transform 



Pd{x) 



exp 



D-l 

V^D{D - 1) ^„2'=/c! ' 



t 



2VD^ 



(S.15) 



In dS.lSb . the expression inside the summation is a 2{D — 1) degree polynomial of t. The power <" is replaced by the Hermite 



polynomial of (— l)"iJ„ ( --^=== j . The derivation of (IS. 15b is based on the Laplace transform relation of C [d"//di"] = s^C[f] 

2 2 

and d"e~^ /dx" = (— l)"i/„(x)e~^ , where C denotes the Laplace transform. 

Figure |5] shows the distribution of D x pjj{x), corresponding to the correlation function in ||2], as compared with K£i{x, x) 
without constraint. The ripple in Dp£,{x) is far larger than the corresponding function without zero- trace constraint. With the 
increase of dimension D, the distribution approaches semicircle distribution, similar to the conclusion in U]. 
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In I?], the monotone sequences in random words are found to have the same statistics as zero-trace Gaussian unitary 
ensemble. Using the notation here, pd{x) is given by the integration equatiorQ: 

—Kd{x,x) = \ — I e~^^ pd{x - y)dy (S.16) 

The original integration equation (IS.16t is for the largest eigenvalue but the same argument applies to poix). The left hand 
side of ( IS.16b is the eigenvalue distribution of Gaussian unitary ensemble without constraint. Using Fourier transform to both 
side of ( IS.16b . Pd{x) can be found and should be the same as ( IS.15b . 

In summary, we have extended the p.d.f. for zero-trace Gaussian unitary ensemble to very high order 
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'^This argument is from Prof. Folkmai' Bomemann of Tech. Univ. Munich, see |http://mathoverflow.net/questions/86965/traceless-gue-four-centered-fermions] 



